function [theta, xi] = f_gmm(y, x, z, W, fe)

    nfe    = size(fe, 2);
    xx     = fe' * fe;
    x1     = x(:, 1:end-nfe);
    xtilda = x1 - fe * (xx \ (fe' * x1));
    ytilda = y  - fe * (xx \ (fe' * y));

    XZ     = xtilda' * z;
    YZ     = ytilda' * z;
    theta1 = inv(XZ * W * XZ') * XZ * W * YZ';

    yfe    = y - x1 * theta1;
    theta2 = inv(xx) * fe' * yfe;

    theta  = [theta1; theta2];
    xi     = y - x * theta;
end
